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ABSTRACT 

We present time-dependent, numerical simulations of the magneto-centrifugal model for jet 



> 

^— N ' formation, in an axisymmetric geometry, using a modification of the ZeusSD code adapted to 

Psl . parallel computers. The gas is supposed cold with negligible thermal pressure throughout. The 

CN I number of boundary conditions imposed on the disk surface is that necessary and sufficient to 
take into account information propagating upstream from the fast and Alfven critical surfaces, 

^\ . avoiding over-determination of the flow and unphysical effects, such as numerical "boundary 

(-H I layers" that otherwise isolate the disk from the fiow and produce impulsive accelerations. 
^ i' It is known that open magnetic field lines can either trap or propel the gas, depending upon 

Q . the inclination angle, 9, of the poloidal field to the disk normal. This inclination is free to adjust, 

H I changing from trapping to propelling when 9 is larger than 9c ~ 30°; however, the ejected mass 

^ ' flux is imposed in these simulations as a function of the radius alone. As there is a region, near 

C^ ■ 

the origin, where the inclination of field lines to the axis is too small to drive a centrifugal wind, 

we inject a thin, axial jet, expected to form electromagnetically near black holes in AGN and 

^^ ■ Galactic superluminal sources. 

H I Rapid acceleration and coUimation of the flow is generally observed when the disk field 

configuration is propelling. We parametrize our runs using a magnetic fiux ^P ex Rr'^'* and mass 

flux j = pvz ex i?~'^J . We show in detail the steady state of a reference run with parameters 

e* = —1/2, Cj — 3/2, finding that the wind leaves the computational volume in the axial 

direction with an Alfven number Ma ~ 4, poloidal speed Vp ^ 1.6w_ffo, collimated inside an 

angle 9 ^ 11°. We show also the thrust T, energy L, torque G and mass discharge M of the 

outgoing wind, and illustrate the dependence of these quantities with the exponents e,f and ej. 

Subject headings: ISM: jets and outflows - galaxies: active - methods: numerical - MHD 



1. Introduction 

Astrophysical jets are fast and well-coUimatcd flows, observed in a wide variety of astronomical 



systems, ranging from young stellar objects (e.g. Lada 1985), Galactic superluminal sources (e.g. Mirabel & 



Rodriguez 1998, HjcUming 1997), to quasars (e.g. Begelman, Blandford & Rees 1984). A single mechanism 
may be responsible for their formation and collimation (e.g. [Livio 1997 ). A candidate model is the 
magneto-centrifugal mechanism proposed by Blandford and Payne (1982), investigated for instance in 



Clarke, Norman fc Burns (1986), Lovelace ct al. (1986), Pudritz fc Norman (1986), Konigl (1989), Ostriker 



(1997), and reviewed recently by Konigl & Pudritz (1999; see, however, Fiege & Henriksen 199(; for an 



alternative model). In this picture, parcels of cold gas are stripped from the surface of a Keplerian disk and 
flung out along open magnetic field lines by centrifugal force. At large distances from the source region, 
rotation winds the field lines up into concentric loops around the axis. These magnetic loops pinch on the 
outflow (or wind) and coUimate it into a narrow jet. 

The original wind solution of Blandford & Payne (1982) is axisymmetric, steady-state and self-similar. 
The self-similarity ansatz is expected to break down near the rotation axis and at large distances, where 
appropriate boundary conditions must be imposed. There have been several studies seeking non-self-similar 



steady-state wind solutions. Sakurai (1985, 1987) and Najita fc Shu (1994) took the approach of starting 
directly from the time-independent MHD equations and found steady-state solutions iteratively. However, 
no general solutions are available at the present from this approach, due to the great complexity of the 
critical surfaces whose loci are unknown a priori (Heinemann & Olbert 1978). Moreover, there is no 
guarantee that a steady-state solution can be found at all for a given set of boundary conditions. A 
time-dependent approach provides a surer method for studying MHD outflows, both steady state or not. 



This more flexible approach has been adopted by several authors (e.g., Lind ct al. 1989; Stone & Norman 



1993|; [Ouyed fc Pudritz 19971 ; [Romanova et al. 1997| ; [Ustyugova et al. 1999| ; [Kcppens fc Goedbloed 199£ 



Bogovalov fc Tsinganos 1999 ). It is also the approach we shall take. 



Following Stone & Norman (1993) and Ouyed & Pudritz (1997), we treat the Keplerian disk as a 
boundary and use the Zeus MHD code to simulate the time-dependent, disk-wind. The main subtleties 
come from treating the disk- wind boundary -the disk surface- and the region close to the rotation axis. A 
special treatment of the disk- wind boundary is warranted, because this is the place where matter is ejected 
into the wind and where magnetic field lines are anchored. On this boundary, some of the flow variables 
can be fixed, but not all of them. This point is illustrated most clearly by steady-state MHD winds, which 
have three, well-known critical surfaces — slow, Alfvenic and fast (Heinemann & Olbert 1978). Crossing 



each of these critical surfaces imposes a set of constraints on the flow variables. Similarly, it is impossible to 
assign arbitrary fixed values to all the quantities at the disk boundary: some quantities must be left free to 
adjust according to flow conditions in the sub-fast region. Failure to do so will overdetermine the problem, 
forming numerical boundary layers with local discontinuities and impulsive acceleration of numerical origin 
close to the disk, which may cast doubts about the origin of the jet launched in the simulation. A numerical 
boundary layer may cause the effective boundary conditions governing the wind flow to differ from those 



imposed at the disk, as mentioned in Meier et al. (1997). This loss of control is particularly inconvenient 
when doing parametrical studies. In this paper we prevent such overdeterminations by properly taking into 
account the existence of information going upstream from the critical surfaces. Our main improvement 
over previous ZEUS-based works is to fix only the necessary number of boundary conditions at the disk 
surface. Romanova et al. (1997) and Ustyugova et al. (1999) have implemented a similar method for their 
Godunov-type code independently. 

In addition, we derive a physically motivated method for treating the outffow near the rotation axis. 
Centrifugal acceleration in the axial region does not happen in the magneto-centrifugal picture because the 
centrifugal force vanishes on the axis. A gravitational infall is expected along the axis, which could disrupt 
the outflow from the disk through streaming instabilities. To prevent such infall, some previous authors 



(e.g., Ouyed & Pudritz 1997) assumed a high (thermal and/or turbulent) pressure to counterbalance the 



gravity of the central object. We, on the other hand, assume the presence of a narrow jet near the axis as 



explained in §2.2. Such a jet could be driven either electromagnetically by a rotating black hole (Blandford 



fc Znajek 1977) in the context of active galactic nuclei (AGN) or by interaction of a central star with a disk 



stellar wind (e.g. Shu et al. 1991) in the context of young stellar objects. 



In this paper we describe in detail our numerical approach and present our first results using the 
code. In §y we explain the physical model, the numerical code is introduced in §3J^ and the the boundary 



conditions in §3^. Our results are described in §y, and conclusions in §|5|. Further features of the numerical 



approach are presented in an Appendix. 

In forthcoming papers we shall explore the effects of the size of the jet-launching disk, the dependence 
of the mass flux j with the inclination angle 9, and perform a 3D study. 



2. Physical Model 

2.1. Equations and notation 

We will use axisymmetric cylindrical coordinates (2, R, (j>), such that the launching surface is at z = 
and the axis at i? = 0. A subindex p will denote quantities in the poloidal plane {z,R). The notation is 



similar to that in Blandford and Payne (1982) or Lovelace et al. (1986) 
The adiabatic MHD equations used here are 

|+V.(pv) = 
/3^+p(v-V)v = -Vp + pV$g+j xB/c 



aB 

'dt 



Vx(vxB) = Vx£: 



du 

at 



V • (uv) = — pV • V 
P = (7-1)" 



where 



p = matter density 

V = velocity flow field 

B = magnetic field 

j = (c/47r)V X B = current density 

£ = V X B = -cE 

<^g — gravitational potential 

p — thermal pressure 

u = internal energy density (per unit volume) 

7 = adiabatic index 

We use a smoothed gravity, $g = ^^o'^q/aA's + -^^ + R^, with VLq and Vg fixed parameters. The product 
WA'o = ^oTg gives the scale of Keplerian speeds. 



Some quantities are conserved along the fieldlines in the steady state (Mcstcl 1968), becoming fmictions 
of the magnetic flux ^ alone. Two such quantities are the relation between mass and magnetic flux 
k = AnpVp/Bp and the conserved angular velocity fl = R~^ {v^ — B^Vp/Bp), which at z = has the 
equilibrium value Q = ^J R^^d'^g/dR = vk/ R- In steady state, the velocity is related to the magnetic field 
by V = ^- + _R51(/). Conservation also applies to the specific angular momentum £ = R {v^ — B^/k) and 
the specific energy e = v'^ /2 + h + ^g — il.RB^/k, where h is the enthalpy, of small importance in our cold 
winds. Finally, the entropy S is conserved by our equations, but it is not relevant here. We will check these 
conserved quantities as a measure of the final approach to steady state. 



2.2. Magneto-centrifugal acceleration model 

Consider poloidal magnetic fieldlines emerging from the surface of an accretion disk. In steady state, 
the flow will be along these lines, due to the frozen-in property of ideal MHD flows. We can picture one 



element of gas flowing along a fieldline as a 'bead' sliding on a rigid 'wire' ( Hcnriksen fc Rayburn 1971 ). Let 
that bead start from z = 0, at a slow speed or at rest. The wire is rotating at a constant angular speed fl. 
This will exert on the bead a centrifugal force outwards, competing against the gravitational force pushing 
the bead inwards. For a lever arm large enough, the centrifugal force wins this competition, and the field 
line propels the fiow. This happens when the inclination 6 of the poloidal fieldlines from the axis is larger 
that a critical angle 9c, allowing the fiow to accelerate centrifugally along the fieldline as if shot by a sling. 
This critical angle is 9c = 30° for Newtonian gravity and Keplerian rotation, only slightly changed by our 
smoothed gravity. 

The fieldlines twist in the toroidal direction due to the differential rotation between the corona 
and the disk. This generates a toroidal component -B^, whose magnetic pressure tends to coUimate 
the fiow, decreasing 9 along a fieldline for increasing height z. In addition, the gradient in the poloidal 
field pressure causes focusing, depending on the shape of the flux function. Critical surfaces are the 
Alfvenic {vp = Bp/ y/Airp = VAp) and fast surface {vp — B/ y/Anp = VAt)- The slow magnetosonic surface 
is not relevant here due to our low sound speed. Acceleration is due to the centrifugal and magnetic 
forces, which taken per unit mass and projected along the poloidal fieldlines, are fc = (v'i/R) sm9 and 
/^f = (j X B)Bp/cpBp = -(l/87r/9i?2)(Bp/Bp)V(i?B0)2 ( [Ustyugova et al. 1999| ). These two forces must 



counter the projected gravitational pull fg = —{'Bp/Bp)V^g. Effective acceleration requires that the 
collimation of the fieldlines is not complete before the Alfven surface is reached: otherwise the projected 
centrifugal force would be too small. 

Close to the axis, the fieldlines have an inclination lower than the critical angle 6c = 30°, so they 
cannot propel the flow centrifugally. For the purpose of designing this fully cold simulation, we cannot 
prevent gravitational infall by the use of thermal pressure, which is an isotropic mechanism that would push 
across fieldlines as well as along them, making a large fraction of our simulation hot. We assume instead 
a ballistic axial injection of matter close to escape speed. Its ram pressure provides a simple non-isotropic 
mechanism, which achieves the numerical objective of stopping gravitational backflow and isolating the 
region where the inclination of fieldlines is too small, without compromising the footpoints of the fieldlines 
with transverse flows, as a pressure-driven mechanism might tend to do. This jet core made of axial flow is 
not only a convenient numerical device; it is physically justified in AGN, and galactic superluminal sources, 
which may create a central electromagnetic jet carrying Poynting flux and exerting an anisotropic Maxwell 



stress at its surface. (Blandford & Znajek 1977). The ballistic injection thus represents the effects of a 



feature of the flow not covered by the magneto-centrifugal model alone. 



3. Numerical Methods 
3.1. Computer code 

To perform our simulations we wrote a variant, called Zeus36, based on the ZeusSD code, released 



by LCA (Clarke, Norman fc Fiedler 1994; Norman 1996; Stone fc Norman 1992). Zeus3D is an explicit, 



finite-differencing algorithm, running on an Eulerian grid. Conservation of mass and momentum components 



is guaranteed by the code. Magnetic fields are evolved using the constrained transport method ( Evans 



fc Hawley 198§| ), which guarantees that no monopoles will be generated, limited only by the numerical 



precision of the machine. An advanced method of characteristics ( Hawley fc Stone 1995 ) is used to calculate 
5 = V X B and the transverse Lorentz force. A large base of tests has already been done on this community 
code, and we have not changed the central difference scheme. 

Zeus36 is a modification of ZeusSD written for parallel computers, such as for instance the Intel 
Paragon and the Touchstone Delta. It uses MPI -a standard message passing interface for parallel machines, 
including networks of workstations. This makes Zeus36 easily portable in the parallel computing world. 
The 3D parallelization scheme is straightforward. We subdivide the computational cubic grid into N smaller 
grids, each given to one of the N compute nodes. Communication between the nodes is realized by taking 
advantage of the boundary condition subroutines, appropriately modified to allow communication between 
adjacent compute nodes. Zeus36 also improves upon the latest released versions of Zeus3D (version 
3.4) by treating the boundary conditions at the corners and edges of the grid more consistently. Corners 
and edges are of less importance when the code is not run in parallel. They become more important in a 
parallel treatment, where they appear in the middle of the overall mesh rather than less obtrusively at the 
boundaries of the calculation. The user interface was modified to make it more convenient to our usage on 
many different computer systems, serial and parallel. For instance, it was decided not to use the powerful 
precompiler Editor ( Clarke 1992| ), bundled into the Zeus3D release. These changes, together with the use 



of dynamic memory allocation, eliminate the need to recompile the main code when the grid size is changed. 

Zeus36 has been tested for both internal consistency and consistency with Zeus3D. We have made 
sure that the output does not depend on the distribution or number of the nodes used to compute a 
given grid. To test consistency with Zeus3D, we have run some test problems using both programs, and 
compared their output byte-by-byte. The small departures observed were those due to the improvement in 
the treatment of corners and edges: no difference appeared after those were accounted for. 



3.2. Boundary conditions 

In these axisymmetric simulations the computational grid has four boundaries: axis [R — 0), outer 
radius (i? = i?max), disk (z = 0) and outer height (z = Zmax)- The code enforces boundary conditions 
by assigning values to the fields p, u, v, and £ at the grid edge and at a few ghost zones outside the 
computational grid. The boundary conditions for B arc enforced indirectly from the values of £, using the 
equation d^/dt = V x £ to evolve B in both active and ghost zones, ensuring V • B = everywhere. The 
exact location and grid staggering of these ghost zones follows Stone & Norman (1992). 



3.2.1. R = 

The axis boundary is geometrically well-determined by axisymmetry, and it has offered no complications 
either in implementing or running. The conditions are implemented as usual for the axis in cylindrical 
coordinates. In the calculation of ghost zone values for R < 0, scalar fields are reflected, and vector 
components are reflected with a change of sign in the directions R and (j). 



O.Z.^. Z — ^max 

The upper boundary is treated in our simulation using an outflow boundary condition, as defined in the 
ZeusSD code. Specifically, for z > Zmax A{z) = A(zniax) for all fields, except for Vz{z) = max(uz(zinax), 0). 
This is not perfect, but, as the physics of our simulation implies the existence of a supersonic, super- Alfvenic 
accelerated outflow, these imperfections are not expected to have a large influence in the simulated flow. 



J./C.o. li — -Kjiiax 

This is more delicate, especially at low altitude where the flow is sub- Alfvenic, and can influence the 
upstream flow. In practice we use similar outflow conditions as for z = Zmax- Label the fieldline that 
becomes trans- Alfvenic at i? = i?max, by $ = ^I^i (Fig. m, and it is clear that artificial conditions imposed 
at i? = i?max for ^ > ^1 may induce spurious coUimation. We must therefore explore the sensitivity of our 
results to our treatment of this part of the flow (cf. |4.3| ). 

In one example of an alternative prescription, Romanova et al. (1997) have used a "force- free" 
prescription, jp||Bp, that can be written as Bp • \J{RB^) = 0; a later proposal from the same authors is 



taking Bp • V{RB^) — aBj^B^, with a a constant to be determined (Ustyugova et al. 1999), observing that 
artificial coUimation may occur for inappropriate box sizes. 



3.2.4. 2 = 

Even more important is the treatment of ^ = 0. It represents the launching surface of the model, 
located at the top of the disk, or the base of the corona, above the slow magnetosonic surface. Its properties 
will determine all the flow downstream. 

The first question to address is the number of boundary conditions we are allowed to fix at this 
surface. This number is equal to the number of waves outgoing normally from the boundary, equal to the 
number of characteristics. We calculate it by counting the degrees of freedom of the system, subtracting all 
constraints. We have, in principle, seven degrees of freedom: the density, three components of v and B, 
minus the constraint V • B = 0, and also the internal energy, because our simulation, while dynamically 
cold, keeps track of a small internal energy. Each crossing of a critical surface is also a constraint that will 
remove one degree of freedom. In principle there are three such surfaces, slow, Alfvenic and fast; however, 
the initial Vp at z = is already larger than the slow speed, so only two critical surfaces will be crossed. We 
are left with flve independent waves out of the original seven. Five boundary conditions should be imposed. 



A more detailed deduction of this result can be seen for instance in Bogovalov (1997) 



We choose to flx the five fields p, u, £^, En and Vz at the disk -if launching were sub-slow, the 



presence of the slow critical surface would reduce the number of independent boundary conditions 
to the first four. At z = 0, the boundary conditions for £^ and En are derived from the infinite 
conductivity of the disk material, giving £^ — Q and Sr = RVlBz- The condition for £^ = (Vp x Bp) • (f) 
implies that Vp||Bp, vn/Bn — v^J B^ = Vp/Bp. The condition for £ji = v^B^ — v^B^ implies that 
v^ — Rn + B^Vp/Bp = Rfl + B^Vz/Bz- Implementation is done by assigning values to fields at the 
ghost zones 2: < 0. Density, velocity and internal energy are defined at the ghost zones as straightforward 
functions of R, with u = pc^o/ij — 1) for a small constant Cgo- This fixes the mass flux pvz', if our launching 
were subsonic, this quantity should be determined from the crossing of the slow critical surface. 

The magnetic fleld evolves in time following dB^/dt ~ R^^dRS^/dR, dBji/dt ~ —d£^/dz, and 
dB^/dt = d£ii/dz — d£z/dR. Fixing 5^ = at z = keeps Bz{R) constant in that plane, anchoring the 
fieldlines to the disk surface. The angle 6 — a,icta.n{Bii/Bz) between the fieldlines and the vertical evolves 
in the active zones; a similar evolution must be present in the ghost zones. Otherwise, the fieldlines would 
present sharp kinks at the base of the wind z = 0, associated with large currents and localized forces which 
might mask the launching mechanism. An evolving Bn requires an £,p{z) dependence in the ghost zones; 
knowing that £(/, vanishes at z = 0, we choose £^ to be an odd function of z. The resulting antisymmetric 
time dependence of the field Bz and the flux ^P allows 6 to vary in approximately the same way in both 
active and ghost zones, so that the final inclination of the fieldlines is decided by the flow itself and its 
crossing of the critical surfaces, producing a much better steady state than alternative approaches in which 
the fieldline inclination is fixed. Similarly we make £r{z) symmetric around its known value at z = 0, 
defining £r{z) — 2RVt{R)Bz{R) — £r{—z) for z < 0. (An alternative possibility here would be extrapolating 
the fieldlines inside the ghost zones; however, this is somewhat more difficult to implement.) With this 
choice, one term of the time dependence of B^ is even in z. We make the other term also even, by requiring 
£z{z) = £z{—z), which allows this quantity to vary freely. Finally we need values for vr and u^. We will 
take them from vr = BrVz/Bz and v^j, — RCl + B^Vz/ Bz^ which we implement using the values of Bz and 
f2 at z = 0. Following the variation of these quantities inside the ghost zones is also a possible option. 

In some simulations we find occasions -early during the run, long before steady state is approached- 
in which i;^ < at the first active zone close to the disk. In that case we modify the boundary conditions, 
allowing the disk to absorb the backfiow, preventing the numerical artifacts shown in the Appendix. The 
ghost values at z < for the fields p, vr and w^ are taken temporarily from their values at the first active 
zone with z > 0, w^ is set to zero and the internal energy to u = pc^solil ^ !)• Thus backfiow is absorbed, 
and will eventually disappear. 

In this paper we artificially prescribe j = pvz as a function of R alone, although mass loading is 
effective only for 6 > 30°. However, in the steady state simulations presented here this inequality is satisfied 
for most of the disk outside the innermost region, justifying our choices after the fact. We shall return to 
this in a following paper, in which j will be allowed to depend also on 6. 

The functions of radius are parametrized using a combination of exponents and softening radii: 

piR) = po/^{i + {R/rprr^ 

Bz{R) = Bzo/VOTWn)^ 

*(i?) = <^ MBzortY- 2^4 

{ ^{Bzorl)\og{l + {R/nf) if 66 = 2 

Vz[R) ~ (''^zinnor + ^^outcr) 






V.outcr = {UoR^iRY I ^{^ + {R/rvoYY- 

We have used v = 1 (linear) and v — 2 (quadratic). For large i?, Vz oc jj(e"o+i/2) = ij-^t.^ 
j EE pw^ ex i?-(«"+'=p) EE i?-'=^ and * oc i?-('=''~2) ^ ^-e*. 

This implementation of disk boundary conditions uses some values of the fields in the corona; there 
is upstream propagation of information. The velocity at z = is supersonic but sub-Alfvenic; upstream 
propagation of waves is physically expected. A simulation that omits this effect is incomplete, and our 
results are, manifestly, sensitive to the treatment of the surface conditions. There is in principle a risk of 
numerical instability involved in using this information in the code; fortunately, our runs did not show this 
kind of instability. 



3.3. Initial conditions 

Naturally, we are interested in flows that are independent of our starting conditions. It turns out 
that transients generally decay in a few Alfven crossing times, and do not influence the late final flow. 
Nevertheless, they deserve some detailed explanation. 

Initially Bp is defined from its boundary condition value Bz{R) at z = 0, and the initial position of 
the footpoint (z = 0, i? = Rq) of a generic poloidal fieldline passing through a point (z, R). This footpoint 
i?o is found from the equation R^ — {R — miri{\ + z/r^Y'^) Ro ^ Rnigri = 0, where B^q, ri,, Cf,, r^, eg and 
mi are constant parameters. The lengthscale ri determines the location of the initial footpoint Re of the 
critical fieldline ^c with 9 ~ 9^, eg gives the initial curvature and coUimation of the ficldlincs, and m£ sets 
the maximum slope tan9 allowed for the fieldlines, through the prescription tan0 = egR/^ri + R/mi) valid 
at z = 0. If mi is set to infinity, all angles 9 are allowed, and Rq{z, R) = R/{1 + zjrYY^ ■ 

The initial density at z > is given by p{z,R) ~ "ripP^z — 0, i?), with rip < 1. Initial Vz in the 
corona is defined similarly through Vz{z,R) — r]yVz{z = 0,R). The initial value for the internal energy is 

In all but one of our simulations we have set the initial value of i?^ to zero, and then consistently made 
the initial ^^(z = 0,R) ^ vk = i?il(i?). For simplicity, this initial v^ is chosen as independent of z. Finally, 
the radial component of the velocity, vji, is set to zero in the bulk of the flow, and to its boundary condition 
value Vfi, = Bji[vz/ Bz) for z < Zmin, where Zmin is the value of z at the first zone above the disk. 



4. Results 

4.1. The reference simulation 

We first describe a fiducial or reference simulation which we shall use to describe the various 
dependencies in our results. In this simulation, we reach the steady state solution shown in Fig. n^. The 
simulation box has 256 x 128 active pixels, with < z < 80.0 and < i? < 40.0. Initial and boundary 
conditions used to produce this run are shown in Fig. 0. The boundary density at z = is constant 



and uniform, p = po = I, with initial p = 0.1 po for z > 0. The Keplerian velocity scale is given by 
vkq = ^Q^g = 1, with Tg = \/3- The function Vz{z — 0, i?) is determined by Vzo = 1-7, fvo — 0.1, 



r, 



9' 



2, Bvo — 1, making e^ — 1.5. Initially Vz{z > 0,R) — 10 ^Vz{R), i?^ = and 



■y^ = Rfl{R), independent of z. The parameters determining Bp are Bzo — 4:VKo\^4:TrpQ, ri, — Vg, e^ = 3/2, 
'''e = T^g, ee = 1 and me = yS = tan 60°. The initial sound speed is set to Cs = 0.0002, ensuring that pressure 
is unimportant in this simulation. The adiabatic index is 7 = 5/3. 

The flow is then evolved until a steady state is reached. Once in steady state, the flow is accelerated 
up to 2.77 times the Keplerian speed at the launching point; the acceleration is purely centrifugal and 
magnetic. The maximum Alfven number Vp/vAp = 3.98 is found close to the upper edge, where also the 
maximum fast Alfven number Vp/vAt is found, equal to 1.36. CoUimation of the fieldlines starts even before 
the flow crosses the Alfven surface Vp = VAp, as we can see in Fig. n^. The minimum Mach number is 
140, showing that pressure is dynamically unimportant in this cold simulation. The acceleration along the 
fieldline '5 = ^e passing through the outer edge of the grid is shown in Fig. ^. Important fieldlines are ^1, 
the last fieldline able to make the Alfven transition; ^c, the fieldline making a critical angle 9c with the 
axis; and vI/q, with footpoint at i? = r„i, at the outer edge of the "core" inner injection, thus separating this 
core from the wind region. 

We have followed this simulation for a time equal to 2.6 Keplerian turns 2Tr/fl a.t R = i?max (and, 
equivalently, 94 turns at the critical radius Re, 170 at the smoothing length radius rg). The poloidal 
velocity and the magnetic field change by less than 1.1 x 10^"* in the last 10% of the run. 

The approach to steady state can also be observed in the figures in the integrals of motion fi, £, /c and 
e, which are already functions of ^ alone, with the exception of border effects for the smallest and largest 
values of z, important mostly for the non-centrifugal core region, ^ < ^Pq- We followed these integrals along 
all fieldlines in the wind ^ > ^Qj and found that on each fieldline the integrals depart from their averaged 
value by less than 5% for fl, 6% for £, 7% for k and 7% for e in the worst case. This still overestimates the 
errors; the corresponding standard deviations divided by the mean are less than 0.3% for fi, and 3% for the 
other integrals, further improved to 0.9% if the outermost region ^ > ^1 is excluded. 

The Alfven radius is defined by R\ — l/fl. Our parameter choice is such that angular momentum is 
primarily extracted from the disk by magnetic stress. Therefore we can estimate £ for low values of z by its 
magnetic term alone, and write an approximate lever- arm ratio as Ra/Rq ~ \/-'{B^/Bz)v\^/vzVk- Fig. ^ 
shows the observed value, the exact formula and the estimate; all agree on a value around 3 for all fieldlines, 
already large enough to justify approximating £ by its magnetic term. In the formula of the estimate, the 
boundary conditions determine all the quantities at z = 0, with the exception of the toroidal field B^, which 
is allowed to self-adjust in the simulation; assuming a value for the lever-arm ratio gives an estimate for the 
toroidal field in steady state. 

The collimation shifts the position of fieldline ^c making a critical angle with the disk surface. In 
steady state 0c is found at Re = 1.85r^, while initially it had been Re = 1.63r£. The integrals M(^), L(vl') 
and G(^) are the total fluxes of mass, energy and angular momentum between the axis and a fieldline ^. 
They have very little numerical dependence on the height z of the surface used to perform the integration, 
because the Zeus algorithm respects these conservation laws. The thrust T is defined here as the surface 
integral of the kinetic energy term pv'^/2, calculated a,t z — Zmax- The integrals M and ^ at 2: = are 
completely determined from the boundary parameters, and have been done analytically. Our results are 
shown in Table |^. 

It is important to observe that in this steady state not only the fieldlines and streamlines collimate. 
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but also the density, which is roughly a function of R alone for large z, as shown in Fig. ^, in spite of the 
boundary condition at z = 0, where p is kept constant. This result is in accord with previous asymptotic 



analysis ( Shu et al. 1995 ) 



As a stringent quality check, we made contour plots of the ratio between the projected magnetic 



and centrifugal forces /m and fc, defined in § |2.2| . Centrifugal acceleration dominates close to the disk, 
approximately up to the point where Vp achieves the local escape speed; magnetic acceleration takes over 
afterwards. There is no visible discontinuity associated with the fieldline \1/ = ^i; this shows that the 
influence of the sub-Alfvenic boundary at i? = i?max is not unduly large. However, we observe a triangle 
near the outer corner where f^ < 0, decelerating the flow -the vertices of this triangle are the points 
{z = 60, R = 40), (z = 80, R = 30) and the outer corner. To decide if this is a physical effect or an artifact 



caused by the edges, we enlarged the box in both directions, as shown below in §4.3.2, finding that it is 
indeed only an artifact. 



4.2. Parametric study 

We now explore the sensitivity of our solution to changes in the functions Vz{R), p{R) and Bz{R), 
parametrized by the exponents e^, Cp and Cf,, or Cj = e„ + Cp and e* = e;, — 2. 

First we checked that Cj is indeed a good parameter. We did a variable density simulation with the 
same value of j = pv^ as the standard run, but with Cp = 0.25 instead of zero, and let it run for the same 
physical time. For z > 2, the results of the variable density simulation coincide with the reference run; 
however, steady state is not achieved at the outer corner (Figs. ^ and 0). Despite its slower convergence, it 
represents the same physical steady state as the reference run, with the exception of very small values of z, 
where the different profile of p still has some influence. For a larger value of the density exponent, such as 
Cp = 1, the run fails to reach steady state in the same box as the reference run. However, simulations with 



a box size four times larger (see §4.3.2) show that the Cp — 1 case is nearly identical to the other cases. The 
effects of the box size are discussed in the next section. We thank the referee for suggesting this test. 

Having done this, we next made a few simulations for different values of the parameters. Fig. H shows 
the critical surfaces and fieldlines for these runs. Their shapes are qualitatively similar, with the remarkable 
exception of the outer fieldlines for the flatter j simulation Cj — 0.5, e^ — 1.5; some of these lines start 
pointing inwards instead of outwards. Not all fieldlines are able to propel the fiow centrifugally at z = 0; 
they must spend some initial kinetic energy before acceleration can start. Further tests can tell if this 
unexpected effect is physical or only a numerical artifact. Because the assumption \9\ > 6c is not fulfilled 
everywhere, these tests should take into account that the mass flux should be reduced for those fieldlines, 
as we will do in future work by using a j dependent on 9. 

Values of speed, Alfven number, lever-arm ratio, luminosity and torque are calculated at the point 
(-^max, ^max), and givcu as functions of the exponents e;,, e^ and Cp in Table y. The first four colums are 
the simulations represented in Fig. y. The last column is the variable density simulation shown in Fig. |6|, 
where we can already see that its flux at the outer corner is slightly different from that of the reference run; 
this is due to its insufficient convergence to the steady state in that region of the computational volume. 
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4.3. Testing of code convergence 



We want to test the convergence of the simulations to a steady state solution by changing some 
computational and physical parameters from our reference run. 



4-3.1. Resolution 

We reduced the resolution of the original standard run by one half, using a box with 128 x 64 active 
pixels, < z < 80.0 and < i? < 40.0. Simulation converged to a steady state similar to the reference run, 
with maximum Vp/vAp — 3.74 and maximum Vp/vAt — 1-33. The resolution test is therefore passed. 



4.3.2. Box size: changes in both length and width together 
Here we show a run done in a larger box with 256 x 128 active pixels, < z < 160.0 and < i? < 80.0; 



this is the same space resolution used in §4.3.1, with a box extending in space to twice the height and width. 
Once a steady state was achieved, we checked that the lower quadrant of this solution was compatible 
with the reference run. For instance, in Fig. S we compare the surfaces of constant Alfven number for the 
two runs with good results: we recover the strongly accelerated flow found before. The maximum Alfven 
number in the lower quadrant of this run is 3.78, and the maximum fast number is 1.34, which is close 



to the reference run and even closer to the results of § 4.3.1 . We also want to check the ratio of projected 
forces /m/ fc, especially near the location of the outer corner of the reference run, now the center of 
the computational box. Looking into Fig. |l^, we see that /m//c in this region is smooth, with positive 
magnetic acceleration; therefore we conclude that the deceleration observed in the reference run was only 
an edge effect related to the outer corner. As another indication that our treatment of the outer boundary 
conditions could be improved for /m, this simulation shows a tiny triangle of magnetic deceleration near its 
outer corner. 



4.3.3. Box size: change in width only; a sub-Alfvenic case. 
This run has also the same physical parameters, changing only the shape of the computational box. 



Here we cut the larger box from §4.3.2 in half, obtaining a square box with 128 x 128 active pixels, 
< z < 80.0 and 0<R< 80.0. 

The results show that the Alfvenic surface crosses the outer boundary z — Zmax^ the flow at the 
boundary R = i?max is fully sub-Alfvenic. The last fieldline able to make the Alfven transition, ^ = ^Pi 
crosses now the z — z^ax edge. It marks the separation between two very different kinds of simulated 
flow. Towards the axis, we have the usual magnetocentrifugally accelerated flow in steady state. But 
outside, in the sub-Alfvenic region, we have a complex structure out of steady state, with filaments of 
magneto-centrifugal acceleration approximately parallel to the fieldlines, wide islands where the projected 
magnetic force /m is negative, and sharp variations in the magnetic twist B^jBp. CoUimation is excessive, 
slowing down the acceleration of the flow. The total projected force along fieldlines, /g + /m + fc, is 
negative in a sizable part of the box, instead of only around the axis. 



As we have shown in §4.3.2, this complex flow disappears in a larger and axially elongated box; it 
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is made possible by the finite size, which sets artificial limits on the critical surfaces. The fieldlines with 
VP < Vf*! are able to reach the Alfven surface, carry the proper critical point information, and propagate it 
downwards to the disk. Due to the shape of the box here, most fieldlines have ^ > ^i, with critical points 
at z > Zjnaxj inaccessible to the simulation. Mathematically, the assumptions used to decide the number of 
functions to fix at the disk surface are wrong, because the critical points are not reached. In this region, 
each streamline is almost independent from each other, creating the possibility of filaments. 

We want to check our assumption that the different behavior of our square and non-square boxes 



is due to the location of the Alfven surface. To do that, we moved the Alfven surface used in § 4.3.3 by 
reducing the intensity of the magnetic fields by half without changing the box. It allows the flow to become 
super- Alfvenic inside the square computational box. The flow is smooth, showing again a case where the 
difference between filamentary and smooth runs lays in the relative position of the simulated Alfven surface 
and the computational volume. We have also checked that the conserved quantities are functions of \1'. 



4-3. 4- Large initial B^: dependence on initial conditions in the sub-Alfvenic case 

The setup used here is similar to the reference run. It differs only in the initial values of B^ and v^, set 
by Vcf,{z = 0, i?) = (1 — fB(t>)vK and B^{z = O.R) = fBc/yVRBz/vz, with /^^ a small fraction, equal to 0.01. 
For simplicity, we have kept using initial values of B^ and v^ independent of z; we could have chosen them 
to make the initial Q a function of ^, but the initial values for B^ would still be away from steady state. 

Here we have a large initial toroidal field, which provides excessive collimation of the fieldlines, making 
the angle 9 too small for magnetocentrifugal acceleration. Also, the magnetic force along the fieldlines is 
such that B = {RB^pY acts like a pressure; B grows along the fieldlines, stopping the acceleration of the 
flow. The results present filamentation, low acceleration, and high twist. This flow is sub-Alfvenic in all the 
volume. We observe here that initial conditions can matter in these simulations, provided that the Alfven 
surface is not reached. 



4-.3.5. Results of the convergence tests 

These tests tell that smooth acceleration of the simulated flow is obtained when a majority of the 
fieldlines reach the critical Alfven surface inside the computational volume. If this condition is met, the 
simulation results depend mostly on the disk boundary conditions, with little dependence on the outer 
boundary conditions. The transients are quickly swept up by the advance of the fast accelerating flow. The 
opposite happens when too many fieldlines lay in the \1/ > ^i region: then the runs show filamentation, 
excessive collimation, dependence on initial and outer boundary conditions, and -obviously- a slow 
acceleration. Controlling the position of the critical surfaces inside the computational box can improve 
convergence, which will help in the setup of future simulations. We plan to utilize this idea in a future work, 
in which we will set up initial and boundary conditions such that the rotating disk will be enclosed by an 
Alfven surface, fully nested inside the computational box; this setup should minimize the infiucncc of the 
outer boundaries. As shown in § 4.3.3, the shape and size of the simulation box can sometimes introduce 
artificial nonsteadiness to a sub-Alfvenic flow. Caution must be observed also in steady-state sub-Alfvenic 
flows, due to their possible dependence on the outer boundary. 

The tendencies to filamentation, both of numerical and physical origin, show that the stability and 
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interaction between the various streamlines in this flow are not fully explored yet. This encourages us even 
more to follow up this work with a 3D simulation, for which this work will be the starting condition. The 
parallel code necessary to do such massive computation is already written and tested. 



5. Conclusions 

This paper is the first in a series exploring various aspects of jets launched from accretion disks. Here 
we have presented the code and shown the results of cold, axisymmetric, time-dependent simulations, 
paying special attention to the treatment of the jet-launching surface and the region near the rotation axis. 
From these simulations we conclude that: 

1) As expected, cold, steady jets can be launched smoothly from Kcplerian disks by the magneto- 
centrifugal mechanism. The same mechanism is also effective in coUimating the outflows. CoUimation is 
observed both in the shape of the fieldlines and in the density profiles, which become cylindrical and thus 
jet-like at large distances along the rotation axis, in agreement with asymptotic analysis. 

2) The magneto-centrifugal mechanism for jet production is robust to changes in the conditions at the 
base of the jets, whose consequences are consistent with expectations. 

3) Steady-state jets obtained in earlier studies using an essentially incomplete treatment of the disk 
boundary conditions are qualitatively unchanged when this deficiency is rectified. It remains to be shown, 
however, whether the same is true for non-steady jets. 

4) The structure of the outflows is insensitive to the initial density or flow speed at the injection surface 
individually, as long as their product (i.e., the mass flux) is kept constant. 

5) The size of the simulation box can strongly influence the outcome of a simulation unless most 
field lines pass through the Alfven surface within the box. As a result, simulations that produce mainly 
sub-Alfvenic disk outflows, steady or not, must be treated with caution. This highlights the pressing need 
for box-invariant simulations, which are under way. 

6) ZeusSD can be parallelized so that solutions from different subgrids match smoothly. 

This research was performed in part using the CACR parallel computer system operated by Caltech. 
We acknowledge useful discussions with James Stone, John Hawley, David Meier and Shinji Koide. This 
research was supported by NSF grant AST-95 29170 and NASA grant 5-2837. RB thanks John BahcaU 
for hospitality at the Institute for Advanced Study and the Beverly and Raymond Sackler Foundation for 
support at the Institute of Astronomy. 



A. Two numerical artifacts related to backflow 

We mention here some numerical problems that had to be overcome in developing the code, in the 
spirit of helping developers of similar codes. A simulation can develop backflow through at least two 
physical mechanisms. One is simple gravitational infall, which operates mostly in the axial region, where 
the magneto-centrifugal acceleration fails. A second mechanism is cocoon backflow. Due to the differential 
Keplerian rotation, the inner parts of the flow launch a jet earlier than the outer. The cocoon of this jet 
may produce a temporary backflow. When the backflow impinges into the disk boundary, our boundary 
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conditions shown in §3.2.4 carefully allow the disk to absorb it. This easy prescription avoids the problems 
discussed below. However, in our earliest simulations, our disk boundary was independent of the sign of Vz 
in the first active zone. This created two kinds of trouble: false acceleration, and filamentation. 



A.l. False acceleration 

In this scenario, backflowing material tries to get back towards the disk, where it finds a boundary 
condition imposing a positive speed outwards. Compression increases the density close to z = 0. This 
increased mass is then flung out at the speed v^ given by the boundary condition imposed. No physical 
process is involved: only a badly chosen boundary. This false acceleration is seen more clearly when we 
reduce the axial injection, allowing the gravitational infall to go all the way into the disk; with incorrect 
boundary conditions, it quickly bounces back unphysically. 



A. 2. Non-physical filamentation 

Let's suppose now that, early in the simulation, a portion of the cocoon backflow touches the disk at 
some radius R = Ri. If the disk boundary condition insists on imposing a forwards speed tj^ > to the 
material, we will get an increase in density due to the artificial compression of the coronal material. If the 
ejection speed Vz is not large, a fast bounce-back will be avoided, and we will not get a false acceleration 
like before. Nevertheless, we still get a localized region, close to the disk, where the density is abnormally 
large. This larger density will be associated with large localized i?0, and it will have a large inertia; the 
magneto-centrifugal acceleration will not be enough to lift the material up. The region around i?i will be 
a hole in the acceleration front; the jet will proceed forwards for both smaller and larger radii, but it will 
stagnate here. A succession of such holes would make the simulation look filamentary; again, allowing the 
disk to absorb the backflow avoids this numerical effect. 
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Fig. 1. — Reference simulation. 

The local Alfven number Vp/vAp is shown in grayscale. In thin black, fieldlines and vectors of Vp. In 
medium black, the critical surfaces Vp = v^p, Vp = v^t and the local escape speed line Vp = -i/— 2$g = Wcsc! 
the line Bp ~ B^ is shown in dashes. In thick black, some selected fieldlines: ^i passing through the 
intersection between the Alfven surface and the outer boundary; V&o at the outer edge of the axial injection 
zone, separating the injected core from the main wind; ^c, the fieldlinc whose inclination is critical at z = 0; 
and ^e, passing through the outer edge of the grid. 
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Fig. 2. — Reference simulation. 

Boundary and initial conditions for Vz^ vaz = Bz/\/4:TTp, VAp, vk = R^, Bz and -v/— $g = v^sd y^ at z ~ 

as functions of i?, normalized by the Keplerian speed scale vkvi = ^o^g- 
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Fig. 3. — Reference simulation. 

Velocities along the fieldline ^ = ^e passing through the outer edge of the grid. Plot of Wp, VApi VAt, ^osc, 

Rfl, RqU and v^ vs. z, showing the growth of poloidal speed. 
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Fig. 4. — Lever arm ratio Ra/Ro (full lines), compared to \/£/fl in dashes and \/—{B^/ Bz)v\^/vzVk in 
dots. Both the flux '^ and the footpoint _Ro are used to label the fieldlines. 
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Fig. 5. — Contours and greyscale of density in the reference run, showing colhmation. 



w=w 



^40 r 




^"S,^ 



p=10" 



p=10 



-25 



p=10" 

p=10 
p=10" 



-15 



-0 



H -0.5 



- -1 



- -1.5 



- -2 



- -2.5 



J -3 



22 



Fig. 6. — A variable density simulation defined by Cp — 0.25, Cj — 1.5, and e^ = —0.5, is compared to the 
reference run, which has the same mass flux j = pvz, and a flat density profile at 2; = 0, defined by Cp = 0.0, 
Bj = 1.5, and e^ = —0.5. Fieldlines and contours of M^ in these two runs are shown. 
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Fig. 7. — The variable density simulation defined in Fig. n is again compared to the reference run, this time 
by showing contours of the density of one run divided by that of the other. 
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Fig. 8. — Comparison of four runs. Alfvcn surface in thick lines, fast surface in medium and ficldlincs in 
thin. Fieldlines have been chosen to pass through the same footpoints in the four runs. 
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Fig. 9. — Contour plots of Ma ~ Vp/vAp, for the reference run (solid lines) and a larger and wider box 
(dashes), compared to check convergence. 
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Fig. 10. — Simulation on the larger and wider box. 

In solid, the contour levels Jul jc = 0; Ij 2, 3, and 4, show the ratio of magnetic over centrifugal forces 

along the fieldlines (shown in dots). The /m = level is a short, thick line at the outer corner. 
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Table 1: Mass discharge M, luminosity L, torque G, and thrust T at selected fieldlines labeled by their 
magnetic flux '5 or footpoint Rq. M is in units of poVKof'^, L in pow^o^^i G in povj^^r^, T in pov^^rl, and 
^I* in units of \/4^TpQVKQrl. Thrust was calculated dX z — Zmax- 

'J'O *c *e 2*e 3^e *1 i?n 



'-max 



L 5.5 13.1 18.1 20.5 24.0 25.0 28.5 

G 8.8 31 54 71 109 127 265 

M 3.2 6.0 7.7 8.6 10.2 10.8 14.8 

* 9.5 22.6 35.8 45.2 67.9 79.2 191 

Ro 1.73 3.20 4.77 6.01 9.41 11.35 40 

Ra/Ro 2.94 3.14 3.29 3.38 3.50 3.51 

T 1.9 4.1 5.5 - - - 



Table 2: Dependence of some flow quantities on the exponents e;,, e^, and Cp. Fluxes and thrust are calculated 
at the outer corner of the grid, and adimensionalized as in Table n^. The maxima of Vp/vKo are taken at 



e„ 


1.5 


0.5 


1.25 


0.5 


1.25 


66 


1.5 


1.5 


1.625 


1.25 


1.5 


Cp 


0.0 


0.0 


0.25 


0.0 


0.25 


ej 


1.5 


0.5 


1.5 


0.5 


1.5 


e^ 


-0.5 


-0.5 


-0.375 


-0.75 


-0.5 


* 


35.8 


42.5 


27.1 


70.8 


38.8 


Ro 


4.77 


5.64 


3.91 


7.68 


5.15 


Ra/Ro 


3.29 


2.31 


3.03 


2.61 


3.31 


Max. Ma 


3.98 


5.06 


5.26 


2.66 


3.78 


Max. Vp/vKo 


1.67 


1.33 


1.67 


1.31 


1.54 


M 


7.73 


15.02 


6.89 


19.54 


8.05 


L 


18.1 


23.5 


14.5 


34.3 


18.8 


G 


54 


85 


38 


162 


58 


T 


5.5 


8.9 


4.8 


10.4 


5.2 



